function TransMat = EstTransMat2 (bin_price, bin_lower, bin_upper, Pt, It, foo, P,I)
% Description:
%   This function estimates the price transition probabilites for each 
%   lagged price and interval. They are estimated using the observed 
%   continuous prices as opposed to the discretized prices. The price 
%   transition probabilites are between discretized prices. 
%
%   P = the number of price bins
%   I = the number of intervals
%   bin_price = the mean price of each defined price bin Px1(mean taken from observed price vector)
%   bin_lower = the lower price in each defined price bin Px1
%   bin_upper = the upper price in each defined price bin Px1
% All information is sorted from lowest bin to highest bin
%
% DATA:
%    Pt = price in each period
%    It = interval in each period
%    Pbint = the bin number of each price in each period
% Data is sorted by date and interval, that is they are sorted in the order
% in which they happened


%disp ('start transition matrix estimation 2')
% get the number of observations
T = size(Pt,1);

% look at the next P
Pnext = Pt([ 2:end  1]);

% generate dummies
Dummies = dummyvar(It);
size(Dummies);
Dummies = Dummies(:,2:end);  % drop one dummy, interval 1

% generate the square of the bins
Pt2 = Pt.*Pt;

% create the interactions
intPt = Dummies.*repmat(Pt,1,I-1 );
intPt2 = Dummies.*repmat(Pt2,1,I-1);

% create the matrix
X = [ones(T,1),Pt,Pt2,Dummies];
%X = [ones(T,1),Pt,Pt2,Dummies, intPt, intPt2];


X = X(1:(T-1),:);  % lop off last obs
Pnext = Pnext(1:(T-1),1);



% % estimate beta
 beta = inv(X'*X)*X'*Pnext;
 uhat = Pnext-X*beta;
 sigma = sqrt(sum(uhat.*uhat)/(size(X,1)));
 
% CONDITIONAL PROBABILITIES Pt+1 | Pt, It

    % calculate expected price for each possible discrete price and interval
        
        % stack discrete prices in order to make a matrix of each possible price/interval combination
        Ptvector = repmat(bin_price, I,1);        
        Pt2vector = Ptvector.*Ptvector;
        
        % make interval dummies in order to make a matrix of each possible
        % price/interval combination
        temp = (1:I)';
        Ivector = temp(:,ones(1,P)).';
        Ivector = Ivector(:);
        D = dummyvar(Ivector);
        D = D(:,2:end);  %cut out the dropped interval
       
        % create the matrix of each possible interval and discrete price
        X = [ones(P*I,1), Ptvector, Pt2vector, D];
        
        % with interactions
%         InteractPt = D.*repmat(Ptvector,1,size(D,2));
%         InteractPt2 = D.*repmat(Ptvector,1,size(D,2));
%         X = [ones(P*I,1), Ptvector, Pt2vector, D, InteractPt, InteractPt2];
%         
        % get the predicted price using estimated parameters
        phat = X*beta;
   
% 
% % calculate the transitition matrix with the normal dist
% % may not calculate correctly for non uniform bins
% % Transmat(Pt, Pt+1,I)
% TransMat = zeros(P,P,I);
% % calculate the probability for each
% for i = 1:I
%     for p = 1:P
%         %(1:P) is a row vector of price bins
%         TransMat(p,:,i) = normpdf(bin_price',phat(p+P*(i-1)),sigma)./sum(normpdf(bin_price',phat(p+P*(i-1)),sigma));
%     end
% end 


% calculate the transition matrix with the emprical dist
TransMat = zeros(P,P,I);
% calculate the probability for each using empricial distribution of the errors

%outside the loop for speed

low_mat = repmat(bin_lower', length(uhat), 1); %replicate in row form for as many price bins as there are
up_mat = repmat(bin_upper', length(uhat), 1); %replicate in row form for as many price bins as there are

% slow because of the sums and repmats within the loop. 10 min with 100
% prices
for i = 1:I
    for p = 1:P
        
        emp_dist = uhat+phat(p+P*(i-1)); % empirical distribution of residuals with the predicted price (Phatt+1 | Pt,It) as the mean
        emp_dist_mat = repmat(emp_dist,1,P); % matrix replication
           
        A = emp_dist_mat>low_mat;
        B = emp_dist_mat>=up_mat;
        num_resids_in_bin = sum(A-B);
        total=sum(num_resids_in_bin,2);  % the count should be the same as the number of residuals

        TransMat(p,:,i) = num_resids_in_bin/total;  % the prob of transitioning to each other price bin
    end
end 
%disp ('end transition matrix estimation 2')






% % TESTING
% for i = 1:I
%     for p = 1:P
%         %(1:P is a row vector of price bins
%         yes(i,p) = sum(TransMat(p,:,i));
%     end
% end
% 
% sumtoone=min(min(yes))



% from matlab help, residuals have different variances which depend on
% the value of their predictors (How?) (heteroskedasticity??)

